/*****************************************************************************************
**		Fast Multipole Algorithm used for gravitational N-body calculation.
**		This code uses an compressed octree data structure
**		and DPMTA's multipole expansion calculation routins.
*****************************************************************************************/
#include "file.h"
#include "tree.h"

//GLOBAL CONSTANTS //////////////////
int Mtermp;               //number of Multipole Expansion terms
int particleLEAF;		//number of particle in a LEAF 
double Wcriteria=1.0;
vektor rootleftbottom;
vektor rootcell;
/////////////////////////////////////

int main( int argc, char *argv[] )
{
	long N = 0;				//number of particle...

	particlePtr partList = NULL;		//particle class array..

	double starttime=0.0;
	double endtime=0.0;
	double stepsize=0.0;
	clock_t start,end;					//used for estimate running time
	char message[400]=" ";				//print the times etc.  to file.
	float timeLoad , timeUpward , timeDownward , timeForce;

/****************** read from filename mass velocities position ************************/
	char filename[50] = "defaultfilenamedata1.dat";			//default file name

	if ( argc == 4 ) 
	{
		strcpy( filename , argv[1] );
		cout<<"reading particle information from file " <<filename<<endl;
		readfilePart( filename , &N , &partList , &starttime , &endtime , &stepsize);
		Mtermp = atoi( argv[2] );
		particleLEAF = atoi( argv[3] );	
	}
		else if ( argc == 5 ) 
	{
		strcpy( filename , argv[1] );
		cout<<"reading particle information from file " <<filename<<endl;
		readfilePart( filename , &N , &partList , &starttime , &endtime , &stepsize);
		Mtermp = atoi( argv[2] );
		particleLEAF = atoi( argv[3] );	
		Wcriteria = atof( argv[4] );
		cout<<"Wcriteria="<<Wcriteria<<endl;
	}
	else
	{
		cout<<"\nparameter problem..\n"
			<<argv[0]<<" N p-term pLEAF \nN is the number of particles, pLEAF is number of particles in LEAF , p-term is the multipole expansion's p-term\n";
		exit(2);
	}

/************* time ****************************************************************************/
#ifdef showtime    		
	start = clock();
#endif    
/*****************************************************************************************/

	cout <<"Starting a COMPRESS FAST MULTIPOLE algorithm "  << " with Lepfrog integration for a \n" << N
     << "-body system,from time t =" << starttime 
	 <<" to t=" << endtime <<" with time step = " << stepsize
	 << " and p = " << Mtermp << " \nand number of max. particle in a LEAF " << particleLEAF << endl;

/************* time ****************************************************************************/
#ifdef showtime    		
	end = clock();
	cout<<"read from file elapsed "<< (end-start) /(float)CLOCKS_PER_SEC<<" seconds.."<< endl;
#endif    
/*****************************************************************************************/

/*************** Calculate Energy ******************************************************
	double startenergy = 0 , endenergy = 0;
	startenergy = calcEnergy( n , partList );
	cout<<"start energy of system = " << startenergy << endl;
****************************************************************************************/
   
	clock_t onestep_time_start;				//all simulation time..

	/* construct tree with the simulation parameter */
	tree fmmtree( partList , stepsize , N );
	
/****************************************************************************************/

	int counter=0;			//number of total step.

	double dt_out= (endtime - starttime) / 10;	//snapshot output frequency.
	double t_out = starttime + dt_out;			//put snapshot on t_out
	start = clock();

	while ( endtime > starttime  )
	{

		onestep_time_start = clock();

		if ( starttime >= t_out) {
			end = clock();
			cout<<"\ntime = " << starttime<<" & stepsize = " << stepsize <<" takes "<<(end-start)/(float)CLOCKS_PER_SEC<<endl;
			t_out += dt_out;

		}
		
/************* time ****************************************************************************/
#ifdef showtime    		
		clock_t start = clock();
#endif    
/*****************************************************************************************/
		float timeSort;

		fmmtree.loadparticles( &timeSort );

/************* time ****************************************************************************/
#ifdef showtime    		
		clock_t end = clock();
		timeLoad = ((end-start)/(float)CLOCKS_PER_SEC);
		cout<<"\nload particles takes "<< timeLoad <<" seconds"<< endl;
		start = clock();
#endif    
/*****************************************************************************************/

		fmmtree.upwardpass();	//	upward pass .. calculate center of mass for the cubes..

/************* time ****************************************************************************/
#ifdef showtime    		
		end = clock();
		timeUpward=( end-start) / (float)CLOCKS_PER_SEC;
		cout<<"\nupward pass ( without load particles) takes "<< timeUpward <<" seconds"<< endl;
		start = clock();
#endif    
/*****************************************************************************************/

		fmmtree.downwardpass( );		// downward pass ... 

/************* time ****************************************************************************/
#ifdef showtime    		
		end = clock();
		timeDownward= ( end-start) / (float)CLOCKS_PER_SEC;
		cout<<"\ndownward pass ( without load particles with upwardpass) takes "<< timeDownward <<" seconds"<< endl;
		start = clock();
#endif    
/*****************************************************************************************/

		fmmtree.calculateforces( );		//	calculate forces with the tree..

/************* time ****************************************************************************/
#ifdef showtime    		
		end = clock();
		timeForce=( end-start) / (float)CLOCKS_PER_SEC;
		cout<<"\ncalculate forces take ( without load particles.."<< timeForce <<" seconds.."<< endl;
#endif    
/*****************************************************************************************/
		
		fmmtree.emptytree( message );
		
		starttime += stepsize;

		counter++;

/************* time ****************************************************************************/
#ifdef showtime          
   clock_t onestep_time_end=clock();
   cout<<"\nall takes "<< ( onestep_time_end-onestep_time_start) / (float)CLOCKS_PER_SEC<<" seconds"<< endl;
   char timetakes[400]=" ";
   sprintf( timetakes , "compFMM;pLEAF;%d;step;%d;p;%d;timeSort;%f;timeLoad;%f;timeUpward;%f;timeDownward;%f;timeForce;%f;all;%f;"\
		     , particleLEAF, counter , Mtermp, timeSort ,  timeLoad , timeUpward , timeDownward , timeForce , ( onestep_time_end-onestep_time_start) / (float)CLOCKS_PER_SEC );
	strcat( timetakes, message );
   /*cout<<message<<endl;*/
   cout<<"number of step = "<<counter;
#endif    
/*****************************************************************************************/


		////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////	
		//	putsnapshotfile( n , partList );
		//	cout<<"calculated forces from the first pass written to the snapshot file.."<<endl;
		////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////	
		putsnapshotfileForce2( N , partList , timetakes );

	}

/************* time ****************************************************************************
#ifdef showtime          
   clock_t allend=clock();
   cout<<"\nall takes "<< ( allend-allstart) / (float)CLOCKS_PER_SEC<<" seconds"<< endl;
   char timetakes[400]=" ";
   sprintf( timetakes , "FMM;p;%d;pLEAF;%d;all;%f;timeLoad;%f;timeUpward;%f;timeDownward;%f;timeForce;%f"\
		    ,Mtermp, particleLEAF,( allend-allstart) / (float)CLOCKS_PER_SEC , timeLoad , timeUpward , timeDownward , timeForce );
	strcat( timetakes, message );
   //cout<<message<<endl;
   cout<<"number of step = "<<counter;
#endif    
/*****************************************************************************************/

	///////////////////////////////////////////////////////////////////////////////
	// now it write file the forces array..
	//putsnapshotfileForce( N , partList , timetakes );

/*************** Calculate Energy ******************************************************
	endenergy = calcEnergy( n , partList );
	cout<<"end energy of system = " << startenergy << endl<<endl;
	cout<<"energy loss = " << startenergy - endenergy;
****************************************************************************************/
	
	cout<<endl<<"end of the FMM"<<endl;
	return 0;
}

